Maps

## Saving 16 x 9.89 in image
## Simple feature collection with 6 features and 19 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -171.0125 ymin: 57.8875 xmax: -101.6625 ymax: 65.5875
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 20
##   City   coordinates `Unnamed: 0` mean_ndvi mean_surface_water mean_uhi_DayNight
##   <chr>  <chr>              <dbl>     <dbl>              <dbl>             <dbl>
## 1 Broch… [-101.6624…          202     1866.               3                   NA
## 2 Eagle  [-141.1958…          305     2408.               2.5                 NA
## 3 Emmon… [-164.5208…          373     2971.              NA                   NA
## 4 Lavre… [-171.0124…          504      794.               2.33                NA
## 5 Lutse… [-110.7375…           91     1330.               2.75                NA
## 6 Mekor… [-166.1958…          337     5845                2                   NA
## # ℹ 14 more variables: mean_uhi_daytime <dbl>, mean_uhi_Nighttime <dbl>,
## #   Continent <chr>, Class <chr>, Skewness <dbl>, Dip_Statistic <dbl>,
## #   Mean_Height <dbl>, shape <fct>, shape_emoji <chr>, shape_color <chr>,
## #   size_emoji <dbl>, lon <dbl>, lat <dbl>, geometry <POINT [°]>
## Warning: There were 2 warnings in `stopifnot()`.
## The first warning was:
## ℹ In argument: `logMean = case_when(...)`.
## Caused by warning:
## ! NaNs produced
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## `summarise()` has grouped output by 'var1'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'var1'. You can override using the
## `.groups` argument.
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -158.0104 ymin: -45.88905 xmax: 178.023 ymax: 69.41286
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 4
## # Groups:   var1 [2]
##   var1                  Shape     Meanvar                               geometry
##   <chr>                 <chr>       <dbl>                       <MULTIPOINT [°]>
## 1 " Daytime UHI (°C)"   Diamond    0.0254 ((-11.73729 7.961042), (15.37292 10.2…
## 2 " Daytime UHI (°C)"   Hourglass  0.818  ((-41.7813 -22.36964), (-42.17087 -22…
## 3 " Daytime UHI (°C)"   Pyramid    0.803  ((-42.52779 -22.27071), (-43.37084 -2…
## 4 " Nighttime UHI (°C)" Diamond    0.522  ((-11.73729 7.961042), (15.37292 10.2…
## 5 " Nighttime UHI (°C)" Hourglass  0.656  ((-41.7813 -22.36964), (-42.17087 -22…
## 6 " Nighttime UHI (°C)" Pyramid    0.487  ((-42.52779 -22.27071), (-43.37084 -2…

## Saving 16 x 9.89 in image
#Violin plots of variables UHI, NDVI, Surface water
ggplot(data = meanval_log, aes(x=Shape, y=logMean , group=interaction(Shape,var1),  fill=var1)) +
  geom_violin()+
  theme_minimal()+
  labs(title = "Violin plot of UHI, NDVI, Surface water mean values " , 
       subtitle = "Based on city shapes")+   
  xlab("City Shape")+
  ylab("Log Mean ")+
  facet_wrap(. ~ var1 , scales = "free", nrow = 2)+
  theme(legend.position = "none",
        panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_rect(fill = "azure"),
        text = element_text(size = 15, face = "bold") ,
        axis.text.x=element_text(angle=30),
        plot.title = element_text(size = 25, face = "bold"),
        plot.subtitle = element_text(size = 15, face = "bold")
        )  # _3_4_5

ggsave("../data/11_Output_R/ViolinPlotsOfVariablesLog.png")
## Saving 16 x 9.89 in image
library("ggpubr")
#       order = c("ctrl", "trt1", "trt2"), 

m1_sum <- meanval_sum %>% 
  filter( var1=='Both Day and Night UHI (°C)')

m1 <- meanval %>% 
  filter( var1=='Both Day and Night UHI (°C)')
msel <- dplyr::select(m1, Shape, Mean)

ggline(msel, x = "Shape", y = "Mean",
       order = c("Diamond", "Hourglass", "Pyramid"), 
       add= "mean_se", group = 1,
       ylab = "Mean", xlab = "City Shape")+
  geom_line(data = m1_sum , aes(x=Shape, y=Meanvar, group=1))+
  labs(title = "Anova plot of UHI Day and Night value " , 
       subtitle = "Based on city shapes")

ggsave("../data/11_Output_R/Anova.png")
## Saving 16 x 9.89 in image
res.aov <- aov(Mean ~ Shape, data = msel)
# Summary of the analysis
summary(res.aov)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Shape          2     42  21.157   15.82 1.44e-07 ***
## Residuals   3745   5008   1.337                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(Mean ~ Shape, data = msel)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Mean by Shape
## Kruskal-Wallis chi-squared = 35.428, df = 2, p-value = 2.027e-08
pairwise.wilcox.test(msel$Mean, msel$Shape,
                 p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  msel$Mean and msel$Shape 
## 
##           Diamond Hourglass
## Hourglass 0.0047  -        
## Pyramid   0.0060  5e-07    
## 
## P value adjustment method: BH
m2_sum <- meanval_sum %>% 
  filter( var1==' Daytime UHI (°C)')

m2 <- meanval %>% 
  filter( var1==' Daytime UHI (°C)')
msel2 <- dplyr::select(m2, Shape, Mean)

ggline(msel2, x = "Shape", y = "Mean",
       order = c("Diamond", "Hourglass", "Pyramid"), 
       add= "mean_se", group = 1,
       ylab = "Mean", xlab = "City Shape")+
  geom_line(data = m2_sum , aes(x=Shape, y=Meanvar, group=1))+
  labs(title = "Anova plot of UHI Day value " , 
       subtitle = "Based on city shapes")

ggsave("../data/11_Output_R/AnovaDay.png")
## Saving 16 x 9.89 in image
res.aov <- aov(Mean ~ Shape, data = msel2)
# Summary of the analysis
summary(res.aov)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## Shape          2      9   4.714   4.082 0.0169 *
## Residuals   3745   4325   1.155                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(Mean ~ Shape, data = msel2)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Mean by Shape
## Kruskal-Wallis chi-squared = 9.4101, df = 2, p-value = 0.00905
pairwise.wilcox.test(msel$Mean, msel$Shape,
                 p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  msel$Mean and msel$Shape 
## 
##           Diamond Hourglass
## Hourglass 0.0047  -        
## Pyramid   0.0060  5e-07    
## 
## P value adjustment method: BH
m3_sum <- meanval_sum %>% 
  filter( var1==' Nighttime UHI (°C)')

m3 <- meanval %>% 
  filter( var1==' Nighttime UHI (°C)')
msel3 <- dplyr::select(m3, Shape, Mean)

ggline(msel3, x = "Shape", y = "Mean",
       order = c("Diamond", "Hourglass", "Pyramid"), 
       add= "mean_se", group = 1,
       ylab = "Mean", xlab = "City Shape")+
  geom_line(data = m3_sum , aes(x=Shape, y=Meanvar, group=1))+
  labs(title = "Anova plot of UHI Day value " , 
       subtitle = "Based on city shapes")

ggsave("../data/11_Output_R/AnovaNight.png")
## Saving 16 x 9.89 in image
res.aov <- aov(Mean ~ Shape, data = msel2)
# Summary of the analysis
summary(res.aov)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## Shape          2      9   4.714   4.082 0.0169 *
## Residuals   3745   4325   1.155                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(Mean ~ Shape, data = msel2)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Mean by Shape
## Kruskal-Wallis chi-squared = 9.4101, df = 2, p-value = 0.00905
pairwise.wilcox.test(msel$Mean, msel$Shape,
                 p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  msel$Mean and msel$Shape 
## 
##           Diamond Hourglass
## Hourglass 0.0047  -        
## Pyramid   0.0060  5e-07    
## 
## P value adjustment method: BH
library(GGally)
library(dplyr)

corval <- cty_sf  |>
  mutate( 
    Shape = case_when( Class == "Diamond"  ~ 1,
                       Class == "Hourglass"  ~ 2  ,
                       Class == "Pyramid" ~  3,
                       Class == "Inverse Pyramid"  ~  4,    
                       TRUE  ~  0 )
    )
      
corval <- dplyr::select(corval,'mean_uhi_DayNight', 'mean_uhi_daytime', 
                        'mean_uhi_Nighttime', 'mean_ndvi' , 'mean_surface_water',
                        'Shape' )
corval <- corval %>% drop_na() 


ggpairs(  corval, columns=c('mean_uhi_DayNight', 'mean_uhi_daytime', 
                            'mean_uhi_Nighttime', 'mean_ndvi' , 'mean_surface_water', 'Shape'))

#ggpairs(  cty_sf, columns=c('mean_uhi_DayNight', 'Class', 'mean_ndvi' , 'mean_surface_water')
#)

ggsave("../data/11_Output_R/pairsDayNight.png")
## Saving 16 x 9.89 in image
ggpairs(  corval, columns=c('mean_uhi_DayNight', 'Shape'))

#ggpairs(  cty_sf, columns=c('mean_uhi_DayNight', 'Class', 'mean_ndvi' , 'mean_surface_water')
#)

ggsave("../data/11_Output_R/pairs.png")
## Saving 16 x 9.89 in image
corval <- dplyr::select(corval,'mean_uhi_DayNight', 'mean_uhi_daytime', 
                        'mean_uhi_Nighttime', 'mean_ndvi' , 'mean_surface_water',
                        'Shape')
corval <- corval %>% dplyr::select(-geometry)
corval <- corval %>% 
  mutate( Shp = Shape
    )
corval2 <- st_drop_geometry(corval)
#cormsel = subset(cormsel, select = -c(Shape,geometry) )
corval2 <- corval2 %>% dplyr::select(-Shape)  # geometry)


cormsel1 <- msel %>% dplyr::select(-geometry)
#cormsel1 = subset(cormsel1, dplyr::select = c(Shape,Mean) ) 
cormsel1 <- cormsel1 %>% 
  mutate( 
    Shp = case_when( Shape == "Diamond"  ~ 1,
                       Shape == "Hourglass"  ~ 2  ,
                       Shape == "Pyramid" ~  3,
                       Shape == "Inverse Pyramid"  ~  4,    
                       TRUE  ~  0 )
    )

cormsel2 <- st_drop_geometry(cormsel1)
#cormsel = subset(cormsel, select = -c(Shape,geometry) )
cormsel2 <- cormsel2 %>% dplyr::select(-Shape)  # geometry)
cormsel2 = subset(cormsel2, select = c(Shp, Mean) ) 
head(cormsel2)
## # A tibble: 6 × 2
##     Shp  Mean
##   <dbl> <dbl>
## 1     3 -6.50
## 2     3 -5.11
## 3     3 -4.78
## 4     2 -4.25
## 5     2 -4.00
## 6     2 -3.85
#ggpairs(  cormsel, columns=c('Mean', 'Shp') )
#ggsave("../data/11_Output_R/pairsDayNight_Shape.png")
library("corrgram")
#head(corval)
#vars2 <- c("mean_uhi_daytime","mean_uhi_Nighttime","mean_ndvi", "mean_surface_water")
vars2 <- c("Shp","Mean")
corrgram( cormsel2,  main="UHI and City Shape")

ggsave("../data/11_Output_R/Corrgram_uhi_city.png")
## Saving 16 x 9.89 in image
corrgram( corval2,  main="UHI and all other variables")

ggsave("../data/11_Output_R/Corrgram_uhi_all_variables.png")
## Saving 16 x 9.89 in image
## class       : SpatRaster 
## dimensions  : 43, 68, 2  (nrow, ncol, nlyr)
## resolution  : 0.002694934, 0.002694934  (x, y)
## extent      : 9.849982, 10.03324, 56.98437, 57.10025  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : Aalborg_uhi.tif 
## names       : Aalborg_uhi_1, Aalborg_uhi_2 
## min values  :     -1.186516,    -0.2753552 
## max values  :      2.113822,     0.9442680 
## [1] "Aalborg_uhi_1"
## class       : SpatRaster 
## dimensions  : 52, 82, 1  (nrow, ncol, nlyr)
## resolution  : 0.002245788, 0.002245788  (x, y)
## extent      : 9.850027, 10.03418, 56.98238, 57.09917  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : Aalborg_NDVI.tif 
## name        : NDVI_mean 
## min value   : -1306.000 
## max value   :  7139.435 
## [1] "NDVI"
## class       : SpatRaster 
## dimensions  : 43, 68, 1  (nrow, ncol, nlyr)
## resolution  : 0.002694946, 0.002694946  (x, y)
## extent      : 9.850027, 10.03328, 56.98463, 57.10051  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : Aalborg_SWATER.tif 
## name        : waterClass 
## min value   :          1 
## max value   :          3 
## [1] "SWATER"

## Saving 16 x 9.89 in image

## Saving 16 x 9.89 in image

## Saving 16 x 9.89 in image

## Saving 16 x 9.89 in image

## Saving 16 x 9.89 in image
# Not Used


#Boxplots of variables UHI, NDVI, Surface water
ggplot(data = meanval, aes(x=Shape, y=Mean , group=interaction(Shape,var1),  fill=var1)) +
  geom_boxplot( width =0 )+
  stat_boxplot(geom = "errorbar") +
  geom_point(data = meanval_sum , aes(x=Shape, y=Meanvar, size=2))+
  geom_line(data = meanval_sum , aes(x=Shape, y=Meanvar, group=var1))+
  theme_minimal()+
  labs(title = "Boxplot of UHI, NDVI, Surface water global mean values " , 
       subtitle = "Based on city shapes")+ 
  xlab("City Shape")+
  ylab("Mean ")+
  facet_wrap(. ~ var1 , scales = "free", nrow = 2)+
  theme(legend.position = "none",
        panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_rect(fill = "azure"),
        text = element_text(size = 15, face = "bold") ,
        axis.text.x=element_text(angle=30),
        plot.title = element_text(size = 25, face = "bold"),
        plot.subtitle = element_text(size = 15, face = "bold")
        )  # _3_4_5

ggsave("../data/11_Output_R/BoxPlotsErrorbarOfVariables.png")
## Saving 16 x 9.89 in image
#Boxplots of variables UHI, NDVI, Surface water log values
ggplot(data = meanval_log, aes(x=Shape, y=logMean , group=interaction(Shape,var1),  fill=var1)) +
  geom_boxplot( width =0 )+
  stat_boxplot(geom = "errorbar") +
  geom_point(data = meanval_log_sum , aes(x=Shape, y=Meanvar, size=2))+
  geom_line(data = meanval_log_sum , aes(x=Shape, y=Meanvar, group=var1))+
  theme_minimal()+
  labs(title = "Boxplot of UHI, NDVI, Surface water global log mean values " , 
       subtitle = "Based on city shapes")+ 
  xlab("City Shape")+
  ylab("LogMean ")+
  facet_wrap(. ~ var1 , scales = "free", nrow = 2)+
  theme(legend.position = "none",
        panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_rect(fill = "azure"),
        text = element_text(size = 15, face = "bold") ,
        axis.text.x=element_text(angle=30),
        plot.title = element_text(size = 25, face = "bold"),
        plot.subtitle = element_text(size = 15, face = "bold")
        )  # _3_4_5

ggsave("../data/11_Output_R/BoxPlotsErrorbarOfVariables_log.png")
## Saving 16 x 9.89 in image